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Statistical  continuum  mechanics  is  used  to  predict  the  coefficient  of  thermal  expansion  (CTE)  for  solid 
oxide  fuel  cell  glass-ceramic  seal  materials  with  different  morphology  and  crystallinity.  Two-point  cor¬ 
relation  functions  are  utilized  to  represent  the  heterogeneous  microstructure  morphology  and  phase 
distribution.  The  model  uses  two-point  correlation  functions  in  conjunction  with  local  properties  to  pre¬ 
dict  the  effective  CTE.  Prediction  results  are  comparable  to  experimental  CTE  results.  The  advantage  of 
using  the  statistical  continuum  mechanics  model  in  predicting  the  effective  properties  of  anisotropic 
media  is  shown,  using  the  ability  to  take  the  microstructure  into  consideration. 
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1.  Introduction 

Statistical  continuum  mechanics  approaches  were  developed  to 
predict  properties  of  heterogeneous  media.  Adams  et  al.  [1,2]  devel¬ 
oped  a  framework,  called  “microstructure  sensitive  design”  (MSD), 
to  guide  information  flow  from  materials  scientists  to  mechani¬ 
cal  designers.  This  method  has  been  used  to  design  and  optimize 
materials  microstructure  and  performance  [3-6].  Using  explicit 
mathematical  relationships  between  the  microstructure  and  prop¬ 
erties,  we  can  identify  the  proper  processing  methods  to  develop 
a  material  with  certain  desired  properties.  In  our  previous  stud¬ 
ies,  the  mechanical  behavior  of  heterogeneous  materials  has  been 
studied,  such  as  elastic  [7-10]  and  inelastic  behavior  [3,5,11-13]. 

This  methodology  also  finds  application  on  fuel  cell  industry. 
For  planar  solid  oxide  fuel  cells  (SOFCs),  a  hermetic  seal  is  needed 
between  the  fuel  and  air  sides  of  the  cell.  Seals  are  commonly 
made  of  glass-ceramics.  An  important  property  of  these  seals  is 
the  coefficient  of  thermal  expansion  (CTE),  required  to  be  close  to 
neighboring  components  in  the  fuel  cell  to  reduce  thermal  induced 
stress.  Operated  at  high  temperatures  (800-1 000 °C)  in  thermal 
cycles,  a  SOFC  requires  superior  thermal  performance.  The  seal  is 
bonded  to  several  other  cell  components,  making  a  failure  in  the 
seal  difficult  to  repair.  Therefore,  it  is  necessary  to  develop  capabil¬ 
ity  to  predict  thermal  induced  stress  in  fuel  cell  due  to  difference 
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from  thermal  expansion  coefficient.  The  glass-ceramic  used  in  this 
study  is  G18,  developed  by  the  Pacific  Northwest  National  Labo¬ 
ratory,  Richland,  WA  [14-16].  G18  is  a  multi-phase  material  with 
crystalline  needles  formed  in  a  glassy  matrix.  Crystallinity  increases 
during  aging  from  an  initial  value  of  55%  to  approximately  72%. 
One  advantage  of  using  glass-ceramics  as  sealant  is  the  capability 
to  adjust  CTE  and  corresponding  thermal  performance  by  tailoring 
the  crystallinity. 

The  motivation  for  this  work  is  to  present  a  methodology  to 
link  the  microstructure  and  all  its  details  to  properties  (in  this 
case  CTE).  The  only  other  comparable  technique  is  the  use  of  finite 
element  method  (FEM),  which  is  very  detailed  and  sometimes 
expensive  (both  in  terms  of  time  and  effort).  In  this  methodol¬ 
ogy,  the  microstructure  is  reduced  to  a  set  of  probability  functions 
and  then  linked  to  properties  using  a  statistical  homogenization 
technique.  The  reduction  of  the  microstructure  to  the  probability 
functions  is  very  similar  to  the  mesh  generation  in  FEM  except 
that  here,  we  have  a  mathematical  relation  instead  of  the  actual 
microstructure. 

Previous  research  in  predicting  the  properties  of  G18  has 
mainly  focused  on  the  elastic  modulus  and  creep  behavior  [11-13]. 
Although  the  macroscopic  CTEs  of  G18  aged  for  4  h  and  lOOOh  are 
known,  the  mathematical  relationship  between  the  microstructure 
and  CTE  has  not  yet  been  studied  to  take  into  consideration  mor¬ 
phology  [11].  We  have  previously  predicted  the  effective  thermal 
mechanical  properties  using  homogenization  averaging  schemes, 
based  on  volume  fraction  and  local  properties.  These  models 
included  the  Rosen-EIashin  [14,15],  Turner  [16],  Kerner  [17],  and 
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Fig.  1.  SEM  micrograph  of  G1 8,  aged  for  4  h  at  750  °C  after  initial  sintering  at  850  °C. 

Levin  [18]  models.  The  results  agree  with  experimental  data  in 
isotropic  media  without  taking  into  consideration  anisotropy.  Also, 
in  the  previous  study,  samples  of  different  aging  times  were  not 
studied,  but  rather,  only  one  overall  CTE  value  was  used  to  compare. 
The  evolution  of  property  and  microstructure  was  not  considered. 

Two-point  correlation  functions  have  been  used  to  describe  the 
microstructure  and  predict  properties  in  the  past  by  Li  et  al.  [6]. 
The  use  of  two-point  correlation  functions  allows  the  prediction 
to  use  the  exact  microstructure  morphology,  while  other  homog¬ 
enization  models  are  based  solely  on  volume  fraction.  Predictions 
of  the  conductivity  demonstrated  the  dependency  on  microstruc¬ 
ture  morphology.  This  formulation  has  been  extended  to  predict 
the  CTE  of  heterogeneous  media.  In  this  study,  coefficients  of  ther¬ 
mal  expansion  of  G18,  aged  for  4h  and  lOOOh  are  predicted  by 
a  Green’s  function  solution  to  the  homogenization  relations  using 
two-point  correlation  functions. 

2.  Two-point  correlation  functions 

Two-point  correlation  functions  are  used  in  a  statistical 
homogenization  framework  to  characterize  different  types  of 
heterogeneous  media.  Correlation  functions  can  describe  the 
microstructure  morphology,  including  volume  fraction,  clustering, 
and  spatial  relationship  of  the  components  [5,8,9,19].  Vectors  with 
varying  orientation  and  magnitude  are  sampled  onto  the  micro¬ 
graph  randomly.  In  this  system,  there  are  three  phases:  barium 
silicate  as  phase  1,  the  matrix  as  phase  2,  and  the  pores  as  phase  3. 
Hexacelsian  is  included  in  the  amorphous  phase  due  to  the  simi¬ 
larities  in  CTEs  [20]. 

SEM  micrographs  of  G18  are  shown  in  Figs.  1  and  2.  Several 
phases  can  be  seen  in  these  micrographs.  The  brightest  phase  is 
barium  silicate;  the  dark  gray  needles  are  hexacelsian;  the  lighter 
gray  phase  is  the  glass  matrix,  and  the  black  phases  are  pores. 


The  minor  phases,  such  as  monoclinic  celsian  and  other  solid 
solutions  are  not  observed  in  these  micrographs  and  therefore 
ignored.  Two-point  correlation  function  is  a  probability  function, 
Pij,  where  i  is  the  phase  of  the  initial  point  and  j  is  the  phase  at  the 
vector  end  point.  In  this  three  phase  microstructure,  these  opera¬ 
tions  result  in  nine  probability  variables,  satisfying  the  normality 
conditions  shown  below,  where  v\  is  the  volume  fraction  of  phase 
i: 


PU+P2\=V\,  P2\+P22=V2 

(1) 

Pij=Pji 

(2) 

r  ->  o,  P;,(r)  ->  Vi 
r  ->•  oo,  Pu(r)  Vj2 

(3) 

As  a  result,  three  diagonal  correlation  functions  will  be  enough 
to  represent  an  anisotropic  microstructure  statistically.  The  orien¬ 
tation  of  the  vector  becomes  important  in  anisotropic  cases,  as  the 
two-point  correlation  function  changes  with  respect  to  the  orien¬ 
tation  angle. 

3.  Statistical  continuum  mechanism  on  thermal  expansion 
tensor 

By  assuming  that  there  are  three  phases  in  the  microstructure, 
the  pore,  amorphous  glass  and  crystalline  barium  silicate,  we  will 
use  the  CTEs  of  these  three  phases  to  predict  the  overall  CTE. 

In  our  previous  works,  we  have  developed  framework  to  pre¬ 
dict  mechanical  properties  and  electrical  properties  of  composites 
and  polycrystalline  materials  [4,5,7,9,10].  Here,  the  framework  is 
extended  to  predict  the  CTE.  We  assume  that  the  heterogeneous 
media  is  composed  of  n  constituents  with  different  CTE  tensors, 
oi[  (i  =  1. .  .n)  and  volume  fractions  vif  with  i  being  the  phase.  The 
size  of  inhomogeneity  is  assumed  to  be  significantly  larger  than 
the  electron  free  path. 

Effective  CTE  tensor  aeff  in  the  heterogeneous  media  is  then 
defined  by  the  equation: 

(st(x))  =  aeff  (VTM)  (4) 

The  brackets  (■■■)  denote  the  ensemble  average.  The  detailed 
procedure  to  solve  aejj  is  explained  in  the  following.  To  define  the 
relationship  between  the  localized  CTE,  a(x),  and  the  ensemble 
average  of  the  CTE,  a0,  we  introduce  the  polarized  CTE,  a(x),  such 
that: 

a(x)  =  otQ  +  a(x)  (5) 

The  polarized  thermal  expansion  is  the  difference  from  the 
ensemble  average  at  location  x.  We  define  the  polarized  CTE  tensor 
P(x)  as: 


P(x)  =  a(x)VT(x) 


(6) 


Since  the  thermal  expansion  strain  has  zero  divergence,  we 
have: 


V  •  £t(x)  =  0  =  Oi0v  •  VT(x)  +  V  •  P(x) 


(7) 


This  is  a  set  of  partial  differential  equations  that  can  be  solved 
by  using  Green’s  function: 


I*' 


g(x,  x')  = 


47ra0  x-x' 


(8A) 

(8B) 


Fig.  2.  SEM  micrograph  of  G18,  aged  for  lOOOh  at  750 °C  after  initial  sintering  at 
850  °C. 


where  T0  is  the  applied  temperature,  x  and  x!  are  the  head  and  tail 
points  of  a  vector  in  the  microstructure.  To  obtain  the  temperature 
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Fig.  4.  Two-point  correlation  function  results  for  G1 8  aged  for  1 000  h  with  the  vector 
orientation  0°  and  90°. 


Fig.  3.  Two-point  correlation  function  results  for  G18  aged  for  4  h  with  the  vector 
orientation  0°  and  90°. 


gradient  field  vT,  Eq.  (8A)  is  differentiated. 

Vr(x)  =  V7o  +  J  dx'C(x  -  x')P(x')  (9) 

where  vT0  is  the  applied  temperature  gradient.  Here,  we  describe  a 
numerical  routine  to  perform  the  integration  over  the  Green’s  func¬ 
tion  for  an  ensemble  of  aggregates  in  a  heterogeneous  medium.  The 
derivation  of  the  Green’s  function  solution,  G(x  -x'),  can  be  seen  in 
Ref.  [21  ].  Due  to  the  singularity  of  the  integral  when  x  =  x\  a  spher¬ 
ical  region  around  the  singular  point  is  excluded.  Using  integration 
by  parts  and  divergence  theorem,  the  Green’s  function  G(x-x')  is 
expressed  as: 

G(x  -  x')  =  -DS(x  -  x')  +  H(x  -  x')  (10) 


Figs.  3  and  4  show  the  two-point  correlation  functions, Pn,  P22 
and  p33,  in  G1 8  aged  for  4  h  and  1 000  h,  respectively.  Since  there  is 
no  obvious  oscillation  period  in  the  data,  there  is  no  clustering  in 
the  microstructure.  On  the  other  hand,  correlation  functions  change 
little  with  respect  to  the  orientation  of  vectors  sampled,  showing 
that  the  microstructure  is  very  isotropic. 

The  formulation  for  the  CTE  formulation  has  been  applied  to  the 
glass-ceramic  described  earlier,  as  G18.  The  CTE  of  the  pores  was 
considered  as  zero,  while  the  barium  silicate  and  the  matrix  phases 
are  shown  below: 
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X  10-6  K-1 
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where  D  and  H  are  defined  by  the  following: 


D  =  5 — 1 

3a0 


H  = 


1  3  fin -I 

Attqq  r3 


(IOA) 

(IOB) 


Here  I  is  the  second  order  identity  tensor  and  h  is  the  unit  vector 
of  x  -  x!. 

Substitute  the  definition  of  the  polarized  field  back  to  Eq.  (10): 

VT(x)  =  VT0  +  J  dx'G(x  -  x')  *  a(x')VT(x')  (11) 

The  average  field  for  state  h,  with  the  state  being  the  phase  in 
this  case,  can  be  calculated  from  the  above  equation: 


(VT(x))ft  =  VTo  +  J  dx'G(x-x')*(a(VT0,h(x')))hVT0  (12) 

Since  no  assumption  is  used  on  representing  the  statistical 
distribution  of  components  in  the  heterogeneous  media,  the  appli¬ 
cation  of  this  statistical  continuum  model  on  the  prediction  of 
thermal  conductivity  can  cover  a  broad  range  of  materials  systems. 


4.  Results  and  discussion 

Using  the  two-point  correlation  formulation  described  in  Sec¬ 
tion  2,  two-point  correlations  have  been  used  on  Figs.  1  and  2,  with 
phase  1  as  the  barium  silicate,  phase  2  as  the  matrix  and  hexacelsian 
needles,  and  phase  3  as  the  pores. 


These  CTE  properties  are  effective  from  room  temperature  up  to 
300  °C  for  soda-lime  glass,  and  up  to  550  °C  for  barium  silicate.  At 
higher  temperatures,  the  barium  silicate  CTE  lowers  and  the  glass 
phase  softens  [20]. 

Glass-ceramic  aged  for  4  h,  shown  in  Fig.  1 ,  is  composed  of  49.6% 
barium  silicate,  48.6%  matrix,  and  1 .6%  pore.  After  aging  for  1 000  h, 
shown  in  Fig.  2,  the  glass  is  composed  of  56.0%  barium  silicate, 
40.0%  matrix,  and  4.0%  pore.  As  shown  in  Figs.  3  and  4,  the  cor¬ 
relation  functions  approximate  to  the  volume  fraction  and  regress 
to  the  square  of  value  of  volume  fractions.  The  statistical  contin¬ 
uum  mechanics  method  was  applied  to  predict  the  effective  CTE  of 
G18.  Results  from  the  statistical  continuum  mechanics  model  are 
shown  below: 
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To  compare  other  common  homogenization  models,  the  classi¬ 
cal  upper  and  lower  bounds,  similar  to  the  Voigt  and  Reuss  models, 
were  used,  in  addition  to  a  self-consistent  model,  which  is  described 
further  in  Refs.  [21-23].  These  models  are  much  more  simple  and 
are  based  on  microstructure  assumptions,  volume  fractions,  and 
individual  phase  properties.  The  upper  and  lower  bounds  related  to 
two  extreme  microstructures  where  the  components  are  arranged 
in  iso-strain  and  iso-stress,  respectively.  The  simpler  models  allow 
us  to  see  the  effect  of  the  microstructure  in  the  effective  CTE, 
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Fig.  5.  Predicted  CTE  by  statistical  mechanics  and  by  the  self-consistent  model  of 
G18  at  different  aging  times. 


when  comparing  to  the  statistical  continuum  mechanics  model. 
The  effective  CTEs  predicted  by  these  models  are  shown  in  Fig.  5. 

Generally,  the  CTE  increases  with  aging.  However,  the  CTE  of 
the  1000-h-aged  sample  in  the  x-direction  is  lower  than  that  of 
4-h-aged  sample.  This  is  due  to  the  higher  amount  of  anisotropy 
and  porosity  in  the  1000-h-aged  sample  compared  to  the  4-h-aged 
sample.  The  anisotropy  can  be  seen  by  comparing  the  ratio  of  a33  to 
an,  giving  1.00  in  4-h-aged  sample  and  1.02  in  1000-h-aged  sam¬ 
ple.  This  indicates  slightly  more  anisotropy  in  the  lOOOh  sample. 
Correspondingly,  the  CTE  along  z-direction  in  1000-h-aged  sample 
increased  slightly  in  comparing  with  that  of  the  4-h-aged  sample. 
This  may  be  due  to  slight  anisotropy  in  the  micrograph  with  some 
preference  in  the  z-direction.  Overall,  the  microstructure  is  gener¬ 
ally  isotropic,  demonstrated  in  the  correlation  function  and  CTE 
results.  The  self-consistent  and  statistical  continuum  mechanics 
results  fall  within  the  classical  upper  and  lower  bounds.  In  gen¬ 
eral,  the  self-consistent  model  gives  higher  predictions  than  the 
statistical  continuum  mechanics  model.  The  statistical  continuum 
mechanics  model  takes  into  consideration  the  slight  anisotropy  of 
the  1000-h-aged  sample.  It  would  be  more  accurate  to  compare 
experimental  data  with  effective  CTE  predictions  calculated  from 
the  actual  micrograph  of  the  experimental  data  sample.  The  statis¬ 
tical  continuum  mechanics  model  provides  results  that  can  be  used 
to  predict  the  CTE  of  similar  glass-ceramic  materials  used  for  SOFC 
seals. 

Although  the  difference  between  the  models  is  small,  it  is  neces¬ 
sary  for  SOFC  seals  to  be  closest  in  CTE  to  neighboring  components 
to  stress  in  the  seal.  The  two-point  correlation  model  provides  a 
method  to  relate  the  microstructure  to  the  effective  properties. 
Because  the  CTE  is  one  of  the  most  important  qualities  in  an  SOFC,  it 
is  necessary  to  find  a  model  that  would  be  most  accurate.  Therefore, 
it  may  be  more  beneficial  to  use  the  statistical  continuum  mechan¬ 
ics  model,  over  the  more  simple  models,  which  make  assumptions 
about  the  microstructure.  The  results  show  that  two-point  corre¬ 
lation  functions  can  serve  as  a  tool  for  future  development  of  SOFC 
seals. 

The  CTE  model  was  also  applied  to  these  two  set  of  samples  elon¬ 
gated  along  z-axis  Figs.  6  and  7  show  the  micrographs  of  samples 
which  are  elongated  images  of  samples  aged  for  4h  and  lOOOh, 
respectively.  The  anisotropy  in  the  micrographs  is  demonstrated 
in  the  corresponding  correlation  functions.  The  crystalline  nee¬ 
dles,  pores,  and  matrix  phases  are  all  elongated.  The  influence  of 
microstructure  evolution  on  the  thermal  behavior  is  reflected  on 
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Fig.  6.  (A)  Microstructure  A,  micrograph  of  G18  aged  for  4  h,  100%  elongated  in  the 
z-axis  and  (B)  the  correlation  function  for  microstructure  A. 


the  CTEs  simulated  by  statistical  continuum  mechanics: 
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The  large  difference  in  CTE  along  different  directions  is  due 
to  the  anisotropy  that  was  introduced  to  the  microstructure.  The 
anisotropy  is  reflected  in  the  ratio  of  a33  to  an,  increasing  dra¬ 
matically  from  1.22  in  elongated  4-h-aged  sample  to  1.37  in 
elongated  1 000-h-aged  sample.  There  is  an  increase  from  the  origi¬ 
nal  anisotropy  indicators.  The  results  show  that  the  CTE  has  lowered 
in  the  x-direction,  while  increasing  in  the  elongated  z-direction. 
These  examples  show  that  the  statistical  continuum  model  takes 
into  account  the  morphology  of  the  microstructure,  as  well  as  the 
volume  fraction,  as  opposed  to  other  homogenization  models,  such 
as  the  self-consistent,  which  only  considers  the  volume  fraction.  In 
the  future,  this  method  can  be  used  on  other  heterogeneous  materi¬ 
als  to  predict  the  CTE.  This  model  can  also  be  used  for  polycrystalline 
materials,  which  can  take  each  crystal  as  a  different  phase,  with 
anisotropic  individual  properties. 


3850 


J.  Milhans  et  al  /  Journal  of  Power  Sources  196(2011)  3846-3850 


B  06 


♦  pi  1  90° 
p22  90° 
▲  p33  90° 
Xpll  0° 
Xp22  0° 
p33  0° 


r  /  pixels 


Fig.  7.  (A)  Microstructure  B,  micrograph  of  G18  aged  for  1000  h,  100%  elongated  in 
the  z-axis  and  (B)  the  correlation  function  for  microstructure  B. 


5.  Conclusions 

A  statistical  continuum  mechanics  model  was  developed  to 
predict  the  CTE  of  heterogeneous  media  composed  by  isotropic 
phases.  Two-point  correlation  models  were  used  to  describe  the 
microstructure’s  morphology,  volume  fractions,  and  phase  distri¬ 
bution.  This  model  was  applied  to  predict  the  CTE  of  glass-ceramic 


SOFC  seal  G18.  This  model  can  be  used  on  other  heterogeneous 
materials,  including  anisotropic  materials,  to  predict  the  effec¬ 
tive  CTE.  Results  show  the  change  in  effective  CTE  with  increased 
barium  silicate.  Also,  for  anisotropic  microstructures,  thermal 
expansion  was  greater  in  the  elongated  direction.  To  cross  validate 
the  accuracy  of  statistical  continuum  mechanics,  predicted  results 
were  compared  with  experimental  data  and  prediction  results  from 
other  models.  Reliability  of  statistical  continuum  mechanics  was 
demonstrated. 
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